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1. SUMMARY 


During this report period, software for determining the delay time between two 
waveforms was further developed and tested on real data. The main improvement was use of 
the multitaper method of spectral analysis which has been shown to yield better results than 
standard single-taper procedures. This procedure constitutes a significant improvement over 
all previous studies of cross-spectral analysis based on the use of a single taper. An applica¬ 
tion to closely-located Yucca Flat explosions recorded at four stations showed significant 
improvement over single-taper results. Inter-source coherence for Pn was found to be 
significandy greater when the two sources lay along the direction of wave propagation than 
when perpendicular to it. Mean velocities along the four source-receiver paths were found to 
be stable and significantly different, suggesting that propagation velocities may vary consider¬ 
ably from one path to another. It seems therefore that at least a part of the observed location 
errors are due to the assumption of a uniform path-independent velocity model in the compu¬ 
tation of location. Analysis of data from two closely located explosions suggests that the 
cross-spectral method may also be used to determine very small differences (of the order of 
1%) between propagation velocities of regional phases. 
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2. INTRODUCTION 


Location relative to a known reference event is often more useful and precise than abso¬ 
lute event location. The use of relative location compensates for most of the errors arising 
from path effects. Such location requires relative arrival times for pairs of events observed at 
common stations, which can be obtained with high precision for similar events by waveform 
cross-correlation. This estimate of the relative time based on the full waveform is generally 
more precise than that based on rather subjective picks of the initial phase arrival. Conven¬ 
tional cross-correlation slides one waveform past another, and the lag with the optimum corre¬ 
lation coefficient is taken as the relative time. However, the cross-spectral analysis method 
can yield relative delay times with resolution up to an order magnitude better than the seismo¬ 
gram sampling interval. 

Most of the effort during the report period was devoted to developing the software for 
multitaper. An application to the available data from Yucca Flat explosions demonstrated 
significant improvement in relative locations. The cross-spectrum method was also applied to 
regional data from two closely located Pahute Mesa, NTS explosions recorded at two common 
stations. Preliminary results suggest that the analysis can detect very small differences (of the 
order of 1%) in the propagation velocities of regional phases generated by the two shots. A 
knowledge of these differences may be useful in a comparison of the near-source properties of 
closely located explosions. 


l 




3. CROSS-SPECTRAL ANALYSIS WITH MULTITAPER 


Single taper cross-spectral analysis of regional data from Yucca Flat explosions with pre¬ 
cisely known locations provided relative locations that were accurate to within about 1 km. 
However, locations based on arrival times determined by alignment of the first peaks were 
found to be somewhat better (Gupta and Davis, 1992). A possible reason for this disappoint¬ 
ing performance was suspected to be the use of a cosine taper on rather short signal windows. 
We therefore applied the multitaper method of spectral analysis which has been shown to 
yield better results than standard single-taper procedures (Park et al ., 1987; Zhu et al ., 1989). 

We employed the multitaper method, based on a family of tapers which are resistant to 
spectral leakage, by mainly following the procedure described by Park et al. (1987) and Ver¬ 
non et al. (1991). A set of k tapers v k (n,w) are generated which are functions of the time 
series length n and the specified inner bandwidth w. The tapers are the eigenvectors of the 
Toeplitz matrix with associated eigenvalues X k . Following Park et al. (1987), w was selected 
to be 4, and only those tapers with A*. greater than 0.9 were used. With these two restrictions, 
k was found to be equal to 7, independent of the number of points in the signal window. 

Consider two signals A and B represented by the time series f a (t) and f b (t), respectively. 
The signals are windowed (e.g. 64 or 128 samples of data), calibrated, demeaned, and 
detrended, providing the signals X a (n) and X b (n), respectively, where n is the number of sam¬ 
ples in the selected window. Let Y a k (N) and Y b (N) represent the Fourier transforms of 

v k X a (n) and v k X b (n), where N = —f-1, respectively. 

An adaptive spectral estimate of f a (t), S a is formed by 





(1) 


£ld a k (N) Y k (N) I 
_k_ 

E ld a(N )! 2 

k 


where the weights d k , which depend on frequency (N) and on the taper order k, are chosen to 
reduce bias from spectral leakage. An initial spectral estimate is made by assuming 


I Y a 1 2 + I Y a !' 


The first guess for weights is then made by 


d a k (N) = 


'IK s a 


K S a + (1-X k ) X [X a (i)i 2 

i=l 

A new spectral estimate of S a is made by inserting the weights from equation (3) into equa¬ 
tion (1) and using this new value into equation (3) to compute the corresponding new weights. 
This process is repeated until the difference between two successive S a is less than 0.0001. 
Once the weights have been determined, the multitaper cross-spectrum is computed from: 


Y„h = 


X (^r 1 X d k (N) Y k (N) d k (N) Y k *(N) 

k k 


k[£(d a W] 1/z t X ( d bW ] ,/2 

k k 

where * indicates the complex conjugate. If Y ab is expressed as 


the phase spectrum <(>(N) is 


Y ab (N) = Re(N) + i Im(N), 


..... Im(N) 

W N) = tan 

Re(N) 


If the two seismograms are identical in shape but have different amplitudes and are shifted in 
time by x, then 


f b (t) = k f a (t+x), 

and by the application of the shift theorem, the phase spectrum becomes 
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<)>(N) = 2 7t N x (8) 

Thus the delay time x can be obtained simply by fitting a straight line through the phase of 

the cross-spectrum with zero intercept. In fitting this slope, the values are weighted based on 

the coherence y ab : 


Yab = 


IY 


ab 1 


Y aa Y, 


(9) 


aa 1 bb 


The weighting factor, W(N) in the linear fit is based on the Hannon-Thomson processor 
(Knapp and Carter, 1976), i. e. 


W(N) = 


Yab(N) 

1 - Yab(N) 


( 10 ) 


The slope of the line has a continuum of possible values, rather than discrete values, 
with the result that estimates of delay time for highly coherent pairs can actually be an order 
of magnitude more precise than the sample rate. This technique has been applied successfully 
to local sequences of earthquakes with impressively precise locations (Ito, 1985; Fremont and 
Malone, 1987). 

Regional data from six closely spaced underground nuclear explosions at the Nevada 
Test Site (NTS), recorded at the four broadband digital stations, ELK, KNB, LAC, and MNV 
(Table 1 and Figure 1), have been analyzed to determine how the cross-spectral method may 
best be used to improve relative locations. For these explosions, with precisely known loca¬ 
tions, the near-source geological and geophysical properties are also known so that the 
influence of parameters such as shot depth and geological environment can also be investi¬ 
gated. Software for determining the delay time between two waveforms has been developed 
and tested. The relative arrival times between two events recorded at a common station are 
determined with a precision of 0.001 sec. As an example, results from the vertical component 
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Figure 1. Map of six Yucca Flat explosions used in the study. Azimuthal directions and dis 
tances from BASEBALL to the four recording stations are indicated. 







ELK records of the Yucca Flat explosions BASEBALL and PALIZA are shown in Figure 2. 
The digital data are sampled at 42 samples/sec and the signal window is 64 points (about 1.5 
sec). Figure 2 (bottom) shows the input waveforms which may be selected to be of any 
desired duration. The signal window over which a cosine taper applies is also variable. The 
top plot shows the phase of the cross-spectrum varying between - n and + k whereas the 
lower plot shows the coherency, varying between 0 and 1, derived by using a bandwidth of 
2.0 Hz. The delay time of 0.0477 (0.0473) sec is simply the slope of the phase of the cross¬ 
spectrum over the specified range of 0.5-5.0 Hz and is obtained by fitting a slope through 0, 
using a weighting scheme based on coherence of the cross-spectrum (equation 10); the quan¬ 
tity within parentheses denotes one standard deviation. 


TABLE 1 

6 YUCCA FLAT EXPLOSIONS USED IN STUDY 


No. 

DATE 

NAME 

m b 

P-WAVE 

SHOT 

DEPTH 





VEL 

DEPTH 

OF WT 





(KM/S) 

(KM) 

(KM) 

1 

06 Sep 1979 

HEARTS (B) 

5.892 

1.763 

0.6400 

0.5070 

2 

14 Nov 1980 

DAUPHIN (A) 

4.554 

1.420 

0.3200 

0.5800 

3 

15 Jan 1981 

BASEBALL (B) 

5.725 

1.970 

0.5639 

0.5120 

4 

01 Oct 1981 

PALIZA (A) 

5.115 

1.497 

0.4724 

0.5300 

5 

11 Nov 1981 

TILCI (A) 

5.035 

1.600 

0.4450 

0.4940 

6 

12 Nov 1981 

ROUSANNE (B) 

5.458 

1.580 

0.5182 

0.4950 


Note: B and A denote shots below and above the water table, respectively. 


Figure 3 shows results of cross-spectral analysis when the signal window is doubled to 
128 points. A comparison with Figure 2 shows a decrease in the average coherency and an 
increase in the delay time. Results when the cosine tapers in Figures 2 and 3 are replaced by 
multitaper are shown in Figures 4 and 5, respectively. Note the significant improvement in 
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Fit: 0.50 • 5.00 Hz 
Weighted 

Intercept constrained 
Oelay: 0.0477 sec 
Sigma: 0.0473 
Phase unwrapped 



2.00 Hz summation 
Av coherence: 0.930 


FREQUENCY Hz 


842.3 nm 0-P Seis 1 
Start time: 15 Jan 81 


Start Point 

7678 

bz 



ELK BASEBALLelk 
20 : 25: 57.8181 


64 pts 

42.010 samp/sec 
Cosine Taper 
Backup: 10 pts 
Arrival: 

20 :25 : 58.0561 


145.5 nm 0-P Seis 1 ELK PALIZAelk 
Start time: 1 Oct 81 19:0:57.9381 


Corrected Arrival: 
19 : 0 : 58.2239 


Start Point 

7683 

bz 


0.0 5.0 Sec 

PROGRAM CSA v ofs 1 8 Wod Oct 21 15TO 23 1992 

Figure 2. Cross-spectral analysis for measuring delay times based on use of the cosine taper. 
The two input waveforms at ELK (bottom), each 64 points long, provide the phase of the 
cross-spectrum (top) and coherence based on bandwidth of 2.0 Hz for frequencies up to 10 
Hz. Mean slope of the phase spectrum, averaged over the indicated frequency range of 0.5- 
5.0 Hz, provides the delay time used to obtain the corrected arrival time for PALIZA. 
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Figure 3. Similar to Figure 2 for signal windows of 128 points. 
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Backup: 21 pts 
Arrival: 
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Figure 4. Similar to Figure 2 based on use of the multitaper (seven tapers). 
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7 tapers 

Fit: 0.50 - 5.00 Hz 
Weighted 

Intercept constrained 
Delay: 0.0552 sec 
Sigma: 0.0177 
Phase not unwrapped 


FREQUENCY Hz 



Av coherence: 0.669 


FREQUENCY Hz 


1187.8 nmO-P Seis 1 ELK BASEBALlelk 
Start time: 15 Jan 81 20 : 25 : 57.8181 



128 pts 

42.010 samp/sec 
Backup: 21 pts 
Arrival: 
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Sla* t Point 
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392.8 nm 0-P Seis 1 ELK PALIZAelk 
Start time: 1 Oct 81 19:0: 57.9381 



Corrected Arrival: 
19:0:58.4932 


0.0 5 0 Sec 


PROGRAM CSA3 vors I t WodOci2l 14 40 21 1992 


Figure 5. Similar to Figure 2 based on use of the multitaper and signal windows of 128 
points. 
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the determination of delay times as indicated by the improved fits of mean slope and the 
reduced values of standard deviation. 

Assuming BASEBALL to be the reference shot with known location and its first peak at 
each station as the reference time, cross-spectral analysis was carried out to obtain the 
corrected arrival times (such as those indicated on Figures 2, 3, 4, and 5) for each explosion 
recorded at the four stations. Results based on use of the cosine taper and signal windows of 
128 points were first determined. For each explosion, the four arrival times are used as input 
to the LOCATE feature of the Analyst Review Station (ARS) to obtain the epicentral location 
and the origin time. The arrival times (without any delay correction) for BASEBALL 
recorded at the four stations are also used to compute the corresponding ARS location. The 
ARS computed locations and the actual (known) locations of six explosions are shown in Fig¬ 
ure 6. The ARS computed locations for the five shots (paired with BASEBALL) are shifted 
by the amount of shift between the ARS computed and actual locations of BASEBALL. A 
comparison of the shifted locations of each of the five shots with its actual location provides 
the location error associated with each shot (Figure 6). The mean location error is only 1.32 
(0.61) km. Cross spectral analysis based on the use of longer (256 points) windows was also 
carried out for the same six explosions. Using a bandwidth of 1.0 Hz for deriving coherency 
and a frequency range of 0.5-5.0 Hz for computing the delay time, the location results were 
found to be very similar to those in Figure 6 in both the magnitude and azimuthal direction of 
the location errors; the mean location error was 1.22 (0.61) km. 

Locations based on the use of arrival times determined by simple alignment of the first 
peaks at common stations are shown in Figure 7. A comparison with locations derived by 
using the cross-spectral analysis (Figure 6) shows significant improvement, especially for the 
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Figure 6. Location results for six Yucca Flat explosions based on arrival times derived from 
cross-spectral analysis using the cosine taper and signal windows of 128 points. Actual loca¬ 
tions, ARS computed locations, and locations with shift based on BASEBALL are shown; the 
five lines indicate the error between the actual and computed locations. Mean location error 
is 1.318 (0.614) km. 
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37.12 



-116.10 -116.00 

Figure 7. Similar to Figure 6 but derived from simple alignment of the first peaks of all 
explosions recorded at the four recording stations. Mean location error is 1.017 (0.218) km. 


13 









two shots PALIZA and TILCI. The mean location error is only 1.017 (0.218) km. Possible 
reasons for the smaller accuracy of the cross-spectral analysis are the spectral leakage and 
variability for the estimates expected when cosine taper is used on rather short time series 
(Vernon et al ., 1991). 

Location results based on use of the multitaper and signal windows of 64 points are 
shown in Figure 8. There is significant improvement compared to the results with cosine 
taper ( Figure 6); the mean location error is 1.010 (0.221) km. These results appear to be 
impressive if one considers the large epicentral distances (about 200 km to 400 km) and the 
complex geology of the Nevada test site. 

Analysis of regional data from Yucca Flat explosions also provides valuable information 
regarding inter-source coherence of closely located explosions recorded at common stations. 
A comparison of the coherency plots in Figures 2 and 3 or Figures 4 and 5 indicate that the 
shorter signal windows of Pn have significantly better coherency over a much larger range of 
frequency than the longer windows. This implies that the first arrivals in Pn are considerably 
more similar than the later arriving phases, perhaps due to greater contamination by later- 
arriving scattered energy. Furthermore, inter-source coherence for Pn appears to be 
significantly greater when the two sources lie along the direction of wave propagation than 
when perpendicular to it. This is illustrated by comparison of coherence at ELK and KNB for 
two pairs of explosions BASEBALL and ROUSANNE (Figure 9) and BASEBALL and PAL¬ 
IZA (Figure 10). These results were based on use of the multitaper and signal windows of 
128 points. For BASEBALL and ROUSANNE, ELK and KNB are nearly along and perpen¬ 
dicular to the direction of wave propagation, respectively. On the other hand, for BASE¬ 
BALL and PALIZA, ELK and KNB approximate directions perpendicular to and along the 
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Figure 8. Similar to Figure 6 but derived from cross-spectral analysis based on use of the 
multitaper and signal windows of 64 points. Mean location error is 1.010 (0.221) km. 
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Figure 9. Waveform and coherency plots for BASEBALL and ROUS ANNE recorded at ELK 
and KNB showing the inter-source coherency for Pn (over the frequency range of 0.1-3.0 Hz) 
to be larger when the two sources lie nearly along the direction of wave propagation than 
when approximately perpendicular to it. 
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Figure 10. Waveform and coherency plots for BASEBALL and PALIZA recorded at ELK 
and KNB, again showing the inter-source coherency for Pn (over the frequency range of 0.1- 
3.0 Hz) to be larger when the two sources lie nearly along the direction of wave propagation 
than when approximately perpendicular to it. 
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wave propagation. These results, based on use of the multitaper and signal windows of 128 
points, are similar to the inter-sensor coherence of regional phases generally observed to be 
greater along the radial direction than along the transverse direction. An attempt will be made 
to exploit these coherency differences to improve locations. Analysis of the Lg phase showed 
the inter-source coherence to be significantly smaller than for Pn, again similar to the 
observed inter-sensor coherence results for Lg and Pn. Also, coherency appeared to vary 
significantly from one station to another. 

In order to understand what causes the discrepancy between the actual locations and 
those computed from observations of Pn from the Yucca Flat explosions, we examined the 
lateral variations in the Pn velocity under the test site. Pn arrival times from closely spaced 
explosions can be used to estimate average velocities to various stations if the spatial locations 
of the explosions are precisely known and data from at least two shots are available. 

Consider a number of closely located sources, S; recorded at a number of well-separated 
receivers, Rj. Let the average P-wave velocities from the limited source region to the rela¬ 
tively distant receivers be ctj. If ry represents the distance from source S, to the receiver Rj 
and Zj are the origin times of the various sources, then the P-wave arrival times, Ty are given 
by 


Tj: = — + Zj 

J a, 


( 11 ) 


Having data from two or more sources at two or more receivers, we can solve for ctj by using 
the method of least squares. 


Using Pn arrival times, based on alignment of the first peaks, from six explosions (Table 
1) recorded at the four stations, the mean velocities to the four stations ELK, KNB, LAC, and 
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MNV were determined to be 7.051 (0.006), 6.636 (0.007), 6.767 (0.007), and 6.469 (0.008), 
respectively. The quantities within parentheses, denoting one standard deviation of the mean, 
are remarkably small so that the velocity differences are significant and outside the margin of 
error. The large differences in mean velocities (maximum about 9%) along the four source- 
receiver paths suggest that propagation velocities may vary considerably from one path to 
another. Moreover, our results indicate velocities along the north-south direction to be some¬ 
what larger than those along the east-west direction, in general agreement with the tectonics 
of the region and the velocity measurements of Harm sen (1992). It seems therefore that at 
least a part of the observed location errors are due to the assumption of a uniform path- 
independent velocity model. 
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4. COMPARISON OF DATA FROM 1 WO CLOSELY LOCATED SHOTS 


Cross-spectrum analysis has been used to provide a high-precision method for measuring 
temporal velocity changes in the earth’s crust (Poupinet et al ., 1987). The method compares 
the scattered wave portions, or coda, of two highly similar waveforms, generated by a pair of 
nearly identical sources (called doublets), recorded at the same receiver at different times. By 
analyzing progressive relative phase delays between the two doublet signals as a function of 
elapsed time, velocity changes on the order of 0.01% could be measured. The high precision 
of the doublet method is possible because coda waves from colocated sources will sample the 
same regions and generally follow paths that increase in length with elapsed time. Thus a 
small decrease in velocity occurring within the sampled region during the time period between 
the doublet origin times will produce a relative delay that increases with elapsed time over 
some portion of the scattered waves. If the velocity change is pervasive, the relative delay 
may increase linearly over the entire coda. This trend of increasing delay is easier to detect 
and measure than an isolated delay based on a single arrival. Using active, repeatable sources 
instead of natural earthquakes, Roberts et al. (1992) used the doublet method for measuring 
small velocity and attenuation changes in solids. 

The Pahute Mesa explosions BUTEO (12 May 1965, shot depth 696 m, yield 0.7 kt) and 
DURYEA (14 April 1966, shot depth 544 m, yield 65 kt) were closely located (probably in 
the same hole) and generated similar records at two common LRSM stations, KN-UT and 
MN-NV at epicentral distances of 321 and 201 km, respectively (Blandford, 1976). Using 20 
samples/sec data at KN-UT, time delays were measured for overlapping 128-point windows 
incremented by 1 sec. A plot of the mean slopes of the phase spectra of each window, aver¬ 
aged over the frequency range of 0.5-2.0 Hz, versus the elapsed time is shown in Figure 11a. 
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Figure 11. Comparison of the LRSM station KN-UT records of DURYEA and BUTEO 
showing (a) delay time, (b) coherency, and (c) differential bandpassed (0.5-2.0 Hz) log-rms 
amplitude versus time starting about 2 sec before the onset of Pn. Results are based on over¬ 
lapping windows of length 6.4 sec for (a) and (b) and 10 sec for (c) with incremental shift of 
1.0 sec for (a) and (b) and 0.5 sec for (c). Bandpassed waveforms of the two explosions are 
shown at bottom. 
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The corresponding coherency is plotted in Figure lib. In Figure 11a, the mean slope, aver¬ 
aged over the first 50 sec of data, appears to suggest a medium velocity differential of about 
0.6%. Bandpassed (0.5-2.0 Hz) log-rms amplitudes, based on overlapping windows of length 
6.4 sec and incremental shift of 0.5 sec, were also obtained for the two shots (see Gupta and 
Wagner, 1992 for more details). Figure 11c shows the differential (DURYEA - BUTEO) 
log-rms amplitude versus time starting about 2 sec before the onset of Pn. Bandpassed 
waveforms of the two explosions are shown at bottom. 

Results of similar processing of data from the same two explosions recorded at MN-NV 
are shown in Figure 12. The mean slope in Figure 12a, averaged over the first 25 sec of data, 
suggests a medium velocity differential of about 1.0%. These results, although preliminary 
and subject to revision in view of the coherency declining rapidly with time, suggest that the 
cross-spectral analysis can be used to determine very small differences between similar 
waveforms of two closely spaced explosions. 
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5. CONCLUSIONS 


A comparison of relative locations based on application of a single taper versus the mul¬ 
titaper to cross-spectral analysis of regional data from Yucca Flat explosions shows significant 
improvement when the multitaper is used. An accuracy of about 1 km has been achieved; an 
impressive result if one considers the large epicentral distances (about 240 and 410 km) and 
the complex geology around the Nevada Test Site. Inter-source coherence for Pn is found to 
be significandy greater when the two sources lie along the direction of wave propagation than 
when perpendicular to it. Mean velocities along paths from the Yucca Flat source region to 
four stations are found to be stable and significantly different, suggesting that propagation 
velocities may vary considerably from one path to another. This means that at least a part of 
the observed location errors are due to the assumption of a uniform path-independent velocity 
model. Analysis of regional data from two closely located shots suggests that the cross- 
spectral method may also be used to determine very small differences (of the order of 1%) 
between propagation velocities of regional phases from pairs of closely spaced explosions. 
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